In silico design of a multi-epitope vaccine against HPV16/18

Background Cervical cancer is the fourth most common cancer affecting women and is caused by human Papillomavirus (HPV) infections that are sexually transmitted. There are currently commercially available prophylactic vaccines that have been shown to protect vaccinated individuals against HPV infections, however, these vaccines have no therapeutic effects for those who are previously infected with the virus. The current study’s aim was to use immunoinformatics to develop a multi-epitope vaccine with therapeutic potential against cervical cancer. Results In this study, T-cell epitopes from E5 and E7 proteins of HPV16/18 were predicted. These epitopes were evaluated and chosen based on their antigenicity, allergenicity, toxicity, and induction of IFN-γ production (only in helper T lymphocytes). Then, the selected epitopes were sequentially linked by appropriate linkers. In addition, a C-terminal fragment of Mycobacterium tuberculosis heat shock protein 70 (HSP70) was used as an adjuvant for the vaccine construct. The physicochemical parameters of the vaccine construct were acceptable. Furthermore, the vaccine was soluble, highly antigenic, and non-allergenic. The vaccine’s 3D model was predicted, and the structural improvement after refinement was confirmed using the Ramachandran plot and ProSA-web. The vaccine’s B-cell epitopes were predicted. Molecular docking analysis showed that the vaccine's refined 3D model had a strong interaction with the Toll-like receptor 4. The structural stability of the vaccine construct was confirmed by molecular dynamics simulation. Codon adaptation was performed in order to achieve efficient vaccine expression in Escherichia coli strain K12 (E. coli). Subsequently, in silico cloning of the multi-epitope vaccine was conducted into pET-28a ( +) expression vector. Conclusions According to the results of bioinformatics analyses, the multi-epitope vaccine is structurally stable, as well as a non-allergic and non-toxic antigen. However, in vitro and in vivo studies are needed to validate the vaccine’s efficacy and safety. If satisfactory results are obtained from in vitro and in vivo studies, the vaccine designed in this study may be effective as a therapeutic vaccine against cervical cancer. Supplementary Information The online version contains supplementary material available at 10.1186/s12859-022-04784-x.


Background
Cervical cancer is the fourth most common cancer affecting women and is caused by human Papillomavirus (HPV) infections that are sexually transmitted [1]. HPV infection occurs in the epithelium through microscopic wounds making basal cells exposed to the virus [2]. HPV is a member of the Papillomaviridae family. HPVs are classified into five distinct major genera, known as alpha, beta, gamma, mu, and nu [3]. Approximately 200 types of HPV have been identified, which are classified into two high-risk and low-risk groups based on their carcinogenic properties. High-risk types are 16, 18, 31, 33, 35, 39, 45, 51, 52, 56, 58, 68, and 59. It is well known that HPV16 and 18 as the most virulent types throughout the world are among the high-risk genotypes that account for 70% of all cervical cancer cases [4]. Types that are low risk include HPV 6,11,40,42,43,44,53,54,61,72,73, and 81 [5].
HPV genome is a circular double-stranded DNA of about 8000 base pairs that replicates in the host cell nucleus [6]. The HPV genome encompasses three regions, namely non-coding long-control region (LCR), as well as early and late regions. Early genes (E1, E2, E4, E5, E6, E7, and E8) and late genes (L1 and L2) are encoded in the two early and late regions, respectively. There are currently three commercial HPV prophylactic vaccines, Gardasil, Cervarix, and Gardasil 9, that are based on virus-like particles and can prevent high-risk HPV infections [7]. Prophylactic vaccines have been successful in preventing HPV infection and cancer. However, there is a significant population of people at high risk for HPV infection and related diseases worldwide. Therefore, the therapeutic vaccine must be designed and manufactured because prophylactic vaccines are ineffective in eradicating past infections and also do not kill infected cells [8].
The activity of p53, which stimulates the expression of genes involved in cell cycle arrest and apoptosis, is inhibited by E6 via the ubiquitin pathway and with the aid of a cellular protein known as E6-associated protein (E6AP) [9]. E7 oncoprotein inactivates the retinoblastoma (Rb) protein through ubiquitin-dependent degradation [10]. Rb prevents cell cycle progression from G1 to S through binding and inactivating transcription factors, such as members of the E2F family [11]. The majority of HPV tumorigenesis studies focus on the role of E6 and E7 proteins [12]. However, the high-risk HPV E5 protein has been demonstrated to have a significant effect on cellular pathways and signaling in human cell lines [13]. E5 protein was found to be expressed during the early stages of carcinogenesis, and its co-expression with E6 or E7 promotes transformation more than either oncoprotein alone [14]. E5 expression can increase the activity of E6 and E7 and has been shown to cause tumor progression [15]. The E5 can stimulate cancer cell proliferation by forming activating complexes with growth factor receptors, such as the Epidermal Growth Factor Receptor (EGFR), resulting in a proliferative state that lasts for a long time [16,17]. The E5 has a function in decreasing cell death, promoting the accumulation of cells with abnormal DNA genetic mutations, and therefore accelerating the malignancy process [18]. Given all of this, it is strongly recommended that the HPV E5 protein be identified as an HPV oncoprotein and considered as a target for vaccine development [19].
The conventional approach of vaccine development is a time-consuming and expensive process requiring pathogen culture in the laboratory [20]. Reverse vaccinology uses genomic information and compute analysis for the development of vaccines without culturing microorganisms [21]. The current study aimed to use the immunoinformatics approach to develop a multi-epitope vaccine against cervical cancer.

Study design
In this regard, cytotoxic T lymphocyte (CTL) and helper T lymphocyte (HTL) epitopes from the E5 and E7 proteins of HPV16 and 18, as well as the C-terminal fragment of heat shock protein 70 (HSP70) from Mycobacterium tuberculosis were used for vaccine design. Evaluation of the physicochemical properties of the vaccine candidate was performed. Subsequently, the multi-epitope vaccine's secondary and 3D structures were predicted. The vaccine's 3D model was refined and validated after that. In the next step, molecular docking analysis of the refined three-dimensional model of the vaccine construct with toll-like receptor 4 (TLR4) was performed and then the stability of the vaccine was evaluated by molecular dynamics (MD) simulation. Finally, after testing the proposed vaccine's expression efficacy in Escherichia coli (E. coli), the in silico cloning of the construct was conducted. The workflow followed in this study is depicted in Fig. 1. Table 1 shows a list of target proteins for multi-epitope vaccine design, along with their accession numbers.

Prediction and selection of T-cell epitopes
The CTL epitopes are essential for the development of long-lasting cellular immunity and can remove the circulating viruses and cells infected with the virus [22]. The HTL epitopes regulate the adaptive immune responses by inducing T-cell cytokines [23]. In the current study, CTL and HTL epitopes were predicted by NetCTL 1.2 and NetMH-CII 2.2 servers, respectively. We selected 9 and 30 high-score CTL and HTL epitopes, respectively, which were able to bind to more types of MHC alleles. Next, they were evaluated for antigenicity, allergenicity, toxicity, and IFN-γ production induction (only for HTL epitopes). Finally, a total of six CTL (Additional file 1: Table S1) and nine HTL epitopes (Additional file 1: Table S2) were selected for the four target proteins.

Molecular interaction analysis of selected CTL and HTL epitopes with HLA alleles
We used the ClusPro 2.0 server to perform molecular docking analysis to ensure that the selected T cell epitopes interacted with the HLA molecules. The results of molecular docking of CTL and HTL epitopes with HLA-A*02:01 and DRB1*01:01 molecules revealed that among CTL epitopes, the YAWVLVFVY epitope had the lowest energy score (− 1033.7 kcal/mol), while among HTL epitopes, AFTVYVFCFLLPMLL epitope with a docking score of − 1162.8 kcal/mol and FVYIPLFLIHTHARF epitope with a docking score of − 1162.6 kcal/mol had the lowest energy score (Table 2) (Additional file 1: Figs. S1 and S2).

Vaccine construct design
The selected HTL epitopes sequence was checked for overlap with the CTL epitopes sequence. It was found that the FVYIPLFLIHTHARF epitope overlapped the FLIHTHARF epitope. Consequently, the insertion of this CTL epitope in the vaccine construct was ignored. The C-terminal of Mycobacterium tuberculosis HSP70 was used as an adjuvant to achieve maximum immune response. The final structure of the multi-epitope vaccine consisted of 509 amino acids, including one adjuvant, five CTL Fig. 1 Schematic representation of the methodology utilized to construct a multi-epitope vaccine to combat HPV16/18 epitopes, and nine HTL epitopes. The linkers included one EAAAK linker, four AAY linkers, and nine GPGPG linkers (Fig. 2).

Analysis of physicochemical, antigenicity, and allergenicity parameters
The ProtParam server predicted the physicochemical properties of the vaccine construct, which are shown in Table 3. The vaccine's composition included a total of 509   amino acids. Tryptophan was the least (0.2%), while alanine and glycine were the most frequent (11%) residues in the vaccine structure (Fig. 3). The solubility of the proposed vaccine was calculated to be 0.987 by the SOLpro server. The VaxiJen server computed an antigenic score of 0.5425 for the main sequence (without the adjuvant) of the vaccine construct. The final sequence antigenicity (along with the adjuvant) with the bacteria and virus models was estimated to be 0.9426 and 0.5221, respectively. The ANTIGENpro predicted the antigenicity of the main sequence of the vaccine to be 0.509057 and it increased to 0.609880 by attaching the adjuvant sequence.

Secondary structure prediction
The vaccine secondary structure was predicted by the PDBsum. The vaccine construct was composed of 12 helices, 10 helix-helix interacts, 55 beta turns, and 29 gamma turns (Fig. 4).

Prediction, refinement, and validation of the 3D structure
I-TASSER server generated five models for vaccine construct using the threading templates (PDB Hit: 5E85, 2KHO, 1DKX, 1U00, and 1DKZ). The C-scores of models 1-5 were − 3.71, − 2.38, − 2.31, − 2.57, and − 2.59, respectively. The C-scores range from − 5 to 2 for the predicted 3D models, and a model with the highest C-score is a high-confidence model. Therefore, model 3 with a C-score of − 2.31 was selected for the refinement (Fig. 5). The GalaxyRefine server was used to refine the vaccine structure. Model 2 was selected based on various parameters, including GDT-HA (0.9401), RMSD (0.45), MolProbity (2.293), Clash score (12.4), Poor rotamers (1.5), and Rama favored (89.9) (Fig. 5). Zlab and ProSA-web were utilized for the evaluation of the vaccine construct quality. The Ramachandran plot analysis of the unrefined showed 86.375%, 8.759%, and 4.866% of residues in the highly preferred, preferred, and questionable regions, respectively (Fig. 6A). After refinement, the highly preferred, preferred, and questionable regions had 95.134%, 3.65%, and 1.217% of residues, respectively ( Fig. 6B). The overall quality score for protein structures was obtained from ProSA-web. The Z-score of the unrefined model was calculated to be − 7.32 ( Fig. 6C), which reached − 7.36 after refinement (Fig. 6D).

Fig. 3
Composition of amino acids in the vaccine construct. The multi-epitope vaccine consists of 509 residues. Alanine and glycine are the highest and tryptophan is the lowest residues in the vaccine

B-cell epitopes prediction
BCPred server predicted 16 linear B-cell epitopes with a length of 20 amino acids (Table 4), which are illustrated in Fig. 7. The ElliPro server predicted 12    (Table 5).

Molecular docking of the vaccine construct with TLR4
The ClusPro 2.0 server was used to perform molecular docking of the vaccine with TLR4, and this server generated 30 clusters. Cluster 2 had the lowest energy score and was selected as the best cluster (Fig. 9A). The energy score for the cluster was -892.5 kcal/mol. The PDBSum result showed that a total of 40 vaccine residues interacted with 46 TLR4 residues (chain A). There were 35 hydrogen bonds between the residues of TLR4 (chain A) and the vaccine construct (Fig. 9B). Table 6 shows the residues involved in the formation of these hydrogen bonds as well as the lengths of the bonds.   (Fig. 10A). The dynamic behavior of alpha carbon atoms in the structure contains adequate information to study important motions in proteins and reflects the general motions of the structure. Therefore, to investigate motion and structural flexibility, the RMSF of alpha carbon atoms was calculated. The vaccine construct was found to be attached to chains A of TLR4 during the molecular docking step. TLR4 has the same sequence in both chains A and B. To correctly assess the effect of the vaccine binding on chain A flexibility, the RMSF value for both chains was calculated. The flexibility of amino acids 150-450 was the same in chains A and B, while the flexibility at the beginning of chain A (amino acids 25-150) and the terminal part of chain A (amino acids 450 onwards) was lower than chain B. The vaccine's RMSF value was very high in the 210-240 and 450-480 regions (Fig. 10B).

Calculating binding free energies
The MM-GBSA and MM-PBSA methods were used to calculate the binding free energies of TLR4, vaccine, and the TLR4-vaccine complex. The total MM-PBSA binding free energy was − 54,518.71 kcal/mol, − 43,712.38 kcal/mol, and − 98,299.95 kcal/mol for TLR4, vaccine, and the TLR4-vaccine complex, respectively. The total MM-GBSA   Table 7.

Codon optimization and in silico cloning
Back translation and optimization of the codon usage of vaccine sequence in E. coli (strain K12) were conducted by JCat. The GC content and CAI value of the vaccine construct were calculated as 55% and 0.97, respectively. Finally, the DNA sequence of the vaccine was inserted between the two sites of the XhoI and EcoRI restriction enzyme in pET28a ( +) by the SnapGene tool (Fig. 11).

Discussion
Cervical cancer is considered to be a major public health concern [1]. There are now commercial prophylactic vaccines available that have been shown to protect vaccinated individuals against HPV infections. High-risk HPV types can cause precancerous lesions  and malignant transformation in those who haven't been vaccinated or are already infected [24]. To fight existing HPV infections and prevent disease progression in cervical cancer, the development of therapeutic HPV vaccines is crucial. Reverse vaccinology is a novel approach that utilizes computational methods for vaccine design [25]. In recent years, this method has been widely used for designing the novel vaccine candidate against different microorganisms, including Hepatitis B Virus [26], dengue virus [27], human cytomegalovirus [28], Onchocerca volvulus [29], Klebsiella pneumoniae [30], Mycobacterium tuberculosis [31], Helicobacter pylori [32], and Candida auris [33]. Most of the efforts in reverse vaccinology to develop the therapeutic HPV vaccine have been limited to the epitopes prediction of E6 and E7 proteins. Jabbar et al. [34] predicted the E6 and E7 antigenic peptides from HPV16 and 18 using immunoinformatics methods and analyzed their ability to bind MHC-I molecules through molecular docking and MD simulation [34]. Yao et al. [35] predicted 59 and 22 CTL epitopes from E6 and E7, respectively, using the IEDB server for the alleles that were mostly distributed in the world population [35]. In a few studies, only chimeric structures containing the epitopes of E5, E6, and E7 proteins have been designed, and no other bioinformatics evaluations have been performed on the structure. So far, only one investigation has been close to our study, in which, using bioinformatics tools for each of the E5 and E7 oncoproteins of the high-risk HPV types 16, 18, 31, and 45, a chimeric structure was designed separately and their bioinformatics evaluation was completed [36].
In the present study, for the first time, we designed the chimeric structure of the T-cell epitopes of E5 and E7 proteins from HPV16 and 18. The CTL and HTL epitopes  were used as vaccine construction blocks following screening for antigenic, allergenic, and toxicity parameters. To assess the interaction of the selected epitopes with HLA alleles, molecular docking was performed between the selected epitopes and two of the most frequent alleles in the world population, including HLA-A*02:01 and DRB1*01:01 [37]. The interaction of the selected epitopes with HLA molecules indicated that they could play a role in eliciting an immune response. Feltkamp et al. pioneers in HPV epitope studies identified the HPV16-E7 sequence RAHYNIVTF (HPV16-E7 49-57 ) as an MHC-I epitope that can elicit CTL-mediated responses and eradicate established HPV l6-induced tumor cells in mice [38,39], this epitope was also predicted and selected in our study.
The vaccine demonstrates weak immunogenicity if only epitopes are included in the vaccine structure [40]. To resolve this problem, adjuvants should be used along with these vaccines [41]. HSPs are highly conserved molecular chaperons present in both prokaryotes and eukaryotes [42]. HSPs are intracellular protective proteins that contribute to the folding of protein structure and maintenance [43]. The binding of HSP70 to the peptide leads to the further induction of peptide presentation by MHC-I and MHC-II and augments the adaptive immune response [44]. Mycobacterium tuberculosis HSP70 is of particular importance in immunity because it has been shown to activate the innate immune system and elicit specific immune responses against tumors and viral infections [45][46][47][48]. HSP70s have a high degree of homology in total proteins between different species. However, the degree of homology differs between distinct parts of the HSP70 molecule [49]. Using the C-terminal region of HSP70s with a lower degree of homology can prevent potential autoimmune reactions in vaccine design [50]. Based on the contents mentioned above, we attached the C-terminal part of Mycobacterium tuberculosis HSP70 to the vaccine as an adjuvant using an EAAAK linker. Moreover, the selected CTL and HTL epitopes were linked by AAY and GPGPG linkers in the vaccine structure, respectively. The AAY linkers promote epitope presentation and diminish junctional epitopes [51,52]. The GPGPG linkers increase the solubility and allow adjacent domains to function more freely [53]. The vaccine candidate was analyzed for physicochemical properties. The molecular weight of the vaccine construct was calculated to be 54.34 kDa making it an effective target for vaccine development because of a molecular weight of less than 110 kDa and the ease of purification process [54]. The theoretical pI was found to be 5.66, indicating the vaccine's acidic nature. The vaccine's stability index was calculated to be 32.07, indicating that it is stable because it is less than 40 [55]. The aliphatic index was calculated to be 90.12, this high aliphatic index indicates that the vaccine is stable over a range of temperatures [56]. The GRAVY score was − 0.174, the negative value of which reveals the hydrophilic nature of the vaccine allowing it to interact better with water molecules. The vaccine sequence was antigenic both alone and with adjuvant, and the binding of an adjuvant to the N-terminal of vaccine construct improved vaccine antigenicity. In addition, our findings indicated that the proposed vaccine was non-allergenic and soluble.
The three-dimensional structure of the multi-epitope vaccine was generated by the I-TASSER server and refined by the GalaxyRefine server. According to the results of the Ramachandran plot and ProSA-web, the 3D model of the vaccine construct significantly improved after refinement. The Ramachandran plot showed that 86.375% of the residues of the unrefined model were found in the highly preferred region, while this value increased to 95.134% after refining. The z-score of the unrefined model was − 7.32, while this parameter was − 7.36 after refining. The more negative z-score confirms improvement after the refinement of the three-dimensional structure of the vaccine [57].
Several studies have revealed that TLR4 expression is upregulated in cervical cancer cells compared with other TLRs and causes apoptosis resistance [58,59]. Therefore, TLR4 was used to conduct the vaccine's molecular docking study. From the clusters generated at the docking stage, a cluster with the lowest energy (− 892.5 kcal/mol) was selected as the input for MD simulation. The RMSD results showed that the stability of the vaccine-TLR4 complex during the simulation period (40 ns) was satisfactory and TLR4 had reached a stable state earlier than the vaccine in the complex. Evaluation of the RMSF plot related to the vaccine demonstrated that vaccine flexibility is very high in the two regions of 210-240 and 450-480, which could be due to the lack of interaction between the two regions and TLR4. The calculation of the binding free energies revealed that the TLR4-vaccine complex is more stable than its individual components. Codon adaptation was done to enhance the expression of vaccine protein in E. coli (strain K12). The GC content and CAI value are essential parameters for evaluating protein expression levels. The GC content and CAI value of the vaccine sequence were 55% and 0.97, respectively. A GC content of 30-70% is required for higher expression [60], and a CAI value of 0.8-1 is considered good for expression in the target organism [61]. As a result, our findings could be regarded as satisfactory. The results of bioinformatics analyzes of this vaccine candidate were promising. However, further in vitro and in vivo experiments are needed to confirm these findings.

Conclusion
HPV prophylactic vaccines have been successful in preventing HPV infection and cancer, but these vaccines do not have therapeutic effects for people infected with the virus. The vaccine development process is costly and time-consuming. Reverse vaccinology is one of the suggested strategies to reduce the cost and time of the vaccine development process. In this study, we used the reverse vaccinology approach to design a multiepitope vaccine against HPV. The CTL and HTL epitopes of the HPV16/18 E5 and E7 proteins were predicted and then tested for antigenicity, allergenicity, and toxicity. Then the selected epitopes were organized into a chimeric structure along with the appropriate linkers and an adjuvant. The vaccine's physicochemical, solubility, antigenicity, and allergenicity parameters were evaluated, and the vaccine's secondary and three-dimensional structures were predicted. To evaluate the binding affinity, molecular docking was performed between the vaccine's refined structure and the TLR4, followed by MD simulation to confirm the vaccine structure's stability. Finally, the vaccine construct was in silico cloned into the pET28a ( +) after the reverse translation. Although the results of the bioinformatics analysis showed that the proposed vaccine can be effective against HPV, but in vitro and in vivo confirmations of these results are needed. The use of nanoparticle delivery systems may help to increase the efficacy of this vaccine. Sanami et al. BMC Bioinformatics (2022) 23:311 Methods

Retrieving protein sequences
The target proteins' amino acid sequences were obtained in FASTA format from the NCBI database (https:// www. ncbi. nlm. nih. gov/).

Vaccine construct design
The selected CTL and HTL epitopes were utilized to design the vaccine construct. Appropriate linkers ensure that epitopes are properly separated and that amino acid residues are flexible enough to fold into appropriate conformations [35]. AAY linkers were used to connect CTL epitopes, whereas GPGPG linkers were used to connect HTL epitopes. The C-terminal of Mycobacterium tuberculosis HSP70 (amino acids 359-625) as an adjuvant was attached at the N-terminal of the vaccine sequence using the EAAAK linker.

Analysis of physicochemical, antigenicity, and allergenicity parameters
Different physicochemical properties of the vaccine construct, such as amino acid number, molecular weight, theoretical pI, instability index, aliphatic index, and GRAVY were predicted using ExPASy ProtParam (https:// web. expasy. org/ protp aram/) [55]. The average molecular weight of an amino acid multiplied by the total number of amino acids in the protein is a standard method for calculating the molecular weight of a protein. The average molecular weight of an amino acid is 110 daltons (Da) [75]. The isoelectric point is the pH at which a protein has no net charge because the positive and negative charges are equal, making it immobile in a direct current electric field. Knowing a protein's theoretical pI can help you choose and optimize protein purification methods like ionexchange chromatography and isoelectric focusing electrophoresis [76]. The instability index provides an estimate of the stability of proteins in a test tube, and proteins with an instability index of less than 40 are predicted to be stable [77]. The aliphatic index of a protein is defined as the relative volume occupied by aliphatic side chains (alanine, valine, isoleucine, and leucine), and this parameter indicates the protein's heat stability [78]. The GRAVY value is calculated by dividing the sum of all amino acid hydropathy values in a protein sequence by the number of residues in the sequence [79]. The solubility of the proposed vaccine was predicted using the SOLpro server [80]. Furthermore, VaxiJen v2.0 and ANTIGENpro predicted the vaccine's antigenicity. ANTIGENpro (http:// scrat ch. prote omics. ics. uci. edu/) predicts protein antigenicity using the data derived from protein microarray analysis [81]. Allergenicity of the vaccine construct was predicted by AllerTOP v.2.0 server.

Secondary structure prediction
PDBsum (http:// www. ebi. ac. uk/ thorn ton-srv/ datab ases/ cgi-bin/ pdbsum/ GetPa ge. pl? pdbco de= index. html) was used to predict the secondary structure of the proposed vaccine. PDBsum provides information on the protein secondary structure, protein-ligand, protein-DNA interactions, and protein structure quality [82].
Prediction, refinement, and validation of the 3D structure I-TASSER server (https:// zhang group. org//I-TASSER/) predicted the three-dimensional structure of the vaccine construct. I-TASSER predicts 3D models from the amino acid sequence and reports a confidence score (C-score) to assess the predicted models' accuracy [83][84][85]. The GalaxyRefine server (https:// galaxy. seokl ab. org/ cgi-bin/ submit. cgi? type= REFINE) was used to increase the quality of the chosen 3D model. This server is dependent on a refining strategy that has been tested in CASP10. This method first reconstructs and repacks the side chains, and then, relaxes the overall structure through MD simulation [86]. The Zlab server and ProSA-web were used for the validation of the initial and refined models. The Zlab server (https:// zlab. umass med. edu/ bu/ rama/) generates the Ramachandran plot. [87]. The overall quality score for input structure is estimated by ProSA-web (https:// prosa. servi ces. came. sbg. ac. at/ prosa. php). If the calculated score falls outside of a range pretty standard of native proteins, the protein structure is likely to have errors [88,89].

B-cell epitopes prediction
BCPred server (http:// ailab-proje cts1. ist. psu. edu: 8080/ bcpred/ predi ct. html) was used to predict linear B-cell epitopes [90]. In this study, the sequence of the vaccine construct was used as input, and other parameters remained as default. Furthermore, the Ellipro server was used for the prediction of discontinuous B-cell epitopes. The ElliPro server (http:// tools. iedb. org/ ellip ro/) uses residue clustering algorithms along with Tornton's technique to predict discontinuous B-cell epitopes [91].

Molecular docking of the vaccine construct with TLR4
Molecular docking is a computational method used to analyze the binding interactions between the ligand and receptor molecules [92]. ClusPro 2.0 server carried out the docking analysis of the vaccine construct with TLR4 (PDB ID: 4G8A). The refined vaccine model and TLR4 were introduced to the server as ligand and receptor, respectively. The interaction between vaccine residues and TLR4 residues in the docked complex was mapped using the PDBsum database (http:// www. ebi. ac. uk/ thorn ton-srv/ datab ases/ cgibin/ pdbsum/ GetPa ge. pl) [93].

Molecular dynamics simulation
The MD simulation is a valuable model for examining the behavior of ligands and receptors in complexes formed during the docking process [94]. The selected cluster from the molecular docking stage was used as input for MD simulation. The GROMACS 2019.6 software was used to conduct MD simulation [94][95][96] and the ff99SB force field was applied to prepare the input structure. The structure's surface charge was neutralized with sodium and chlorine ions, and protein was introduced into a layer of TIP3P water molecules with a thickness of 10 angstroms using gmx solvate software. To reduce van der Waals interactions and the formation of hydrogen bonds between water and complex molecules, the energy of the structures was minimized using the steepest descent approach. The system temperature was then steadily raised from 0 to 300 K in 200 ps at constant volume, and the system was equilibrated at constant pressure. The MD simulation was conducted at a temperature of 300 K and a period of 40 ns. The calculation of non-bonded interactions with a cutoff of 10 angstroms was performed by the PME algorithms.

Calculating binding free energies
At this stage, the binding free energies for TLR4, vaccine, and the TLR4-vaccine complex were calculated by the MMPBSA method and using the gmx_MMPBSA tool. The calculation of binding free energies was performed using Poisson-Boltzmann (PB) and generalized Born (GB) approaches [97]. The binding energies were calculated throughout the simulation system using 1000 frames selected at regular intervals from simulation trajectories.